Time Series Modeling, Analysis and Decomposition

---

Let me tell you the story of how

I lost $98 Million…

ggplot(corn %>% 
         filter(date > '2007-01-01' & date < '2016-01-01'), 
       aes(x = date, y = value)) + 
  geom_line() + ggtitle('Historical Corn Prices - USD ($) Per Bushel')

Working at a fund

• Futures trading - what’s a future?
• How do you know what a good price is?
  • Build a model, of course… this is where the quants come in
• How do you set a lookback window?
  • Based off of industry, sector, 99% domain knowledge
  • Does it correlate with anything else? – extremely risky!
  • Is the correlation lagged (rain & crops)
  • Knowledge of time series is important
• Setup something to give a confidence interval on prices
  • Make a decision as to how risky you want to be
  • As prices move decide when to buy/sell

ggplot(rain %>%
         filter(date > '2007-01-01'), 
       aes(x = date, y = pr, col=Country)) + geom_line()

c = corn %>% 
  tbl_time(date) %>%
  as_period('monthly', side='start') %>%
  mutate(date = lubridate::floor_date(date,'month'))

r = rain %>%
  select(date, Country, pr) %>%
  spread(Country, pr)

df = c %>%
  left_join(r, by = 'date') %>%
  rename(corn = value) %>%
  drop_na() #%>% gather(nam,val,-date)

p = ggplot(df) 
p + geom_point(aes(x = corn, y = AUS)) 

p + geom_point(aes(x = corn, y = BRA))

p + geom_point(aes(x = corn, y = USA))

plot(df %>% select(-date))

Time series analysis:

• Stationary vs. Non-Stationary
• Additive vs. Multiplicative
• ARMA / ARIMA models
• Filters: Kalman, Holt-Winters
• Stochastic Processes, Random Walk, Wiener Process, Brownian Motion

What is a stationary time series?

“A time series is stationary if a shift in time doesn’t cause a change in the shape of the distribution. Basic properties of the distribution like the mean, variance, and covariance are constant over time.”

Stationary vs Seasonal Time series - http://www.statisticshowto.com/stationarity/

Different types of stationarity:

• Strict Stationarity - basically doesn’t exist… essentially nothing changes based off of time…
• First-Order Stationarity - means never change, but other statistics can change
• Second-Order (Weak) Stationarity - constant mean, variance, and autocovariance that doesn’t change with time
• Trend-Stationary Models - fluctuate around deterministic trend (series mean). Trend is typically linear or quadratic, amplitude of fluctuations stays constant
• Difference-Stationary Models - one or more differences of the data need to be taken to make it stationary

Why do we care about stationarity?

Most statistical models that have been developed assume stationarity.

# Example from: https://www.analyticsvidhya.com/blog/2015/12/complete-tutorial-time-series-modeling/
AirPassengers

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1949 112 118 132 129 121 135 148 148 136 119 104 118
## 1950 115 126 141 135 125 149 170 170 158 133 114 140
## 1951 145 150 178 163 172 178 199 199 184 162 146 166
## 1952 171 180 193 181 183 218 230 242 209 191 172 194
## 1953 196 196 236 235 229 243 264 272 237 211 180 201
## 1954 204 188 235 227 234 264 302 293 259 229 203 229
## 1955 242 233 267 269 270 315 364 347 312 274 237 278
## 1956 284 277 317 313 318 374 413 405 355 306 271 306
## 1957 315 301 356 348 355 422 465 467 404 347 305 336
## 1958 340 318 362 348 363 435 491 505 404 359 310 337
## 1959 360 342 406 396 420 472 548 559 463 407 362 405
## 1960 417 391 419 461 472 535 622 606 508 461 390 432

plot(AirPassengers)
abline(reg = lm(AirPassengers ~ time(AirPassengers)))

#convert to df for rolling operations
ap = matrix(AirPassengers, nco = 1)
t1 = as.integer(min(time(AirPassengers)))
t2 = as.integer(max(time(AirPassengers)))
y = sort(rep(seq(from = t1, to = t2, by = 1), t2-t1 + 1))
m = rep(seq(from = 1, to = 12, by = 1), t2-t1 + 1)
apdf = tibble(year = y, month = m, passengers = ap[,1])
apdf = apdf %>% mutate(date = as.Date(paste0(year,'-',month,'-01'))) %>%
  select(date,passengers)

#remove last 3 years of data
apdf_3 = apdf %>% mutate(x = 1:nrow(apdf)) %>% filter(date < '1958-01-01')
ts_fit = lm(passengers ~ x, data = apdf_3)
ts_fit

## 
## Call:
## lm(formula = passengers ~ x, data = apdf_3)
## 
## Coefficients:
## (Intercept)            x  
##      95.042        2.493

intercept = ts_fit$coefficients[1]
slope = ts_fit$coefficients[2]

apdf_prediction = apdf_3 %>% mutate(type = 'historical') %>%
  bind_rows(apdf %>% filter(date >= '1958-01-01') %>% mutate(type='actual')) %>%
  mutate(x = 1:nrow(apdf)) %>%
  mutate(prediction = (x*slope) + intercept)

ggplot(apdf_prediction, aes(x=date, y=passengers, col=type)) + 
  geom_line() + 
  geom_line(aes(y=prediction), col='black')

rolling_mean = rollify(mean, window = 8)
tmp = apdf %>% mutate(ma = rolling_mean(passengers)) %>% drop_na()
ggplot(tmp, aes(x=date)) + geom_line(aes(y=passengers)) + geom_line(aes(y=ma), col = 'red')

Decomposition Example

Multiplicative model

Forecast = Trend * Seasonality * Cycle * Error

Decomposition Google Sheet

Do transformations!

plot(log(AirPassengers))

plot(diff(AirPassengers))

plot(diff(log(AirPassengers)))

adf.test(diff(log(AirPassengers)))

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(log(AirPassengers))
## Dickey-Fuller = -6.4313, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

acf(log(AirPassengers))

ACF is good for identifying MA component

acf(diff(log(AirPassengers)))

PACF is good for identifying AR component

pacf(diff(log(AirPassengers)))

---